/* SPU Parallel Elliptic Curve Library - Provides Elliptic Curve Cryptographic
 * functions to the Synergistic Processing Unit target platform
 *
 * Copyright (C) 2008  Michael Annichiarico (mike.gcdeveloper@gmail.com)
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
#include"spec.h"
#include<string.h>
#include<libmpm.h>

void nothing() {
	int i=0;
	i+=2;
}

/*
r = 6901746346790563787434755862277025555839812737345013555379383634485463
 == 0x1000000000000000000000000000013E974E72F8A6922031D2603CFE0D7ull
Polynomial Basis:
b = 066 647ede6c 332c7f8c 0923bb58 213b333b 20e9ce42 81fe115f 7d8f90ad
Gx = 0fa c9dfcbac 8313bb21 39f1bb75 5fef65bc 391f8b36 f8f8eb73 71fd558b
Gy = 100 6a08a419 03350678 e58528be bf8a0bef f867a7ca 36716f7e 01f81052
c = sqrt(b) == 00000187f85627b9 7874e747ee31e06d 71caaeea52f21253 e5f946d061da9138
*/
ECGROUP nistb233grp = {
	.n = field2m_init_ull4(10000000000,0000000000000000,0013E974E72F8A69,22031D2603CFE0D7),
	.nbe = field2m_init_ul8be(100,00000000,00000000,00000000,0013E974,E72F8A69,22031D26,03CFE0D7),
	.a = field2m_init_ul8(066,647ede6c,332c7f8c,0923bb58,213b333b,20e9ce42,81fe115f,7d8f90ad),
	.b = field2m_init_ul8(066,647ede6c,332c7f8c,0923bb58,213b333b,20e9ce42,81fe115f,7d8f90ad),
	.c = field2m_init_ull4(00000187f85627b9,7874e747ee31e06d,71caaeea52f21253,e5f946d061da9138),
	.G = {
		.x = field2m_init_ul8(0fa,c9dfcbac,8313bb21,39f1bb75,5fef65bc,391f8b36,f8f8eb73,71fd558b),
		.y = field2m_init_ul8(100,6a08a419,03350678,e58528be,bf8a0bef,f867a7ca,36716f7e,01f81052),
		.z = field2m_init_ull4(0,0,0,1)
	}	
};

void spec_point_print(PECPOINT p) {
	printf("X: ");
	poly_print(&p->x);
	printf("\nY: ");
	poly_print(&p->y);
	printf("\nZ: ");
	poly_print(&p->z);
	printf("\n");
}

ECPOINT spec_point_LDinf = { .x = poly_const_one, .y = poly_const_zero, .z = poly_const_zero };
ECPOINT spec_point_Ainf = { .x = poly_const_zero, .y = poly_const_zero, .z = poly_const_one };

__inline int spec_point_cmpeq(PECPOINT P, PECPOINT Q) {
//	return(memcmp(P,Q,2*sizeof(FIELD2M))==0?1:0);
	return(poly_cmpeq(&P->x,&Q->x) && poly_cmpeq(&P->y,&Q->y)?1:0);
}
__inline int spec_ppoint_cmpeq(PECPOINT P, PECPOINT Q) {
//	return(memcmp(P,Q,3*sizeof(FIELD2M))==0?1:0);
	return(poly_cmpeq(&P->x,&Q->x) && poly_cmpeq(&P->y,&Q->y) && poly_cmpeq(&P->z,&Q->z)?1:0);
}


void spec_point_LD2A(PECGROUP GRP, PECPOINT Q, PECPOINT P) {
	FIELD2M T1;

	if(spec_ppoint_cmpeq(P,&spec_point_LDinf)) {
		*Q=spec_point_Ainf;
		return;
	}
	if(poly_cmpeq(&P->z,&poly_zero)) {
		printf("ERROR LD2A P->z == 0 !!!\n");
		return;
	}
	if(poly_cmpeq(&P->z,&poly_one)) {
		if(Q!=P) *Q=*P;
		return;
	}
	poly_inv_mod(&T1,&P->z);
	poly_mod_mult(&Q->x,&P->x,&T1);
	poly_mod_sqr(&T1,&T1);
	poly_mod_mult(&Q->y,&P->y,&T1);
	Q->z = poly_one;
}

void spec_point_LDdbl(PECGROUP GRP, PECPOINT P3, PECPOINT P1) {
	FIELD2M T1,T2;

	if(spec_ppoint_cmpeq(P1,&spec_point_LDinf)) {
		*P3=spec_point_LDinf;
		return;
	}

	poly_mod_sqr(&T1,&P1->z);
	poly_mod_sqr(&T2,&P1->x);
	poly_mod_mult(&P3->z,&T1,&T2);
	poly_mod_sqr(&P3->x,&T2);
	poly_mod_sqr(&T1,&T1);
	poly_mod_mult(&T2,&T1,&GRP->b);
	poly_add(&P3->x,&P3->x,&T2);
	poly_mod_sqr(&T1,&P1->y);
	poly_add(&T1,&T1,&P3->z);//if a==1 (i think it does)
	poly_add(&T1,&T1,&T2);
	poly_mod_mult(&P3->y,&P3->x,&T1);
	poly_mod_mult(&T1,&T2,&P3->z);
	poly_add(&P3->y,&P3->y,&T1);
}
void spec_point_MLDdbl(PECGROUP GRP, PECPOINT P3, PECPOINT P1) {
	FIELD2M T1,T2;

	if(spec_ppoint_cmpeq(P1,&spec_point_LDinf)) {
		*P3=spec_point_LDinf;
		return;
	}

	poly_mod_sqr(&T1,&P1->z);
	poly_mod_sqr(&T2,&P1->x);
	poly_mod_mult(&P3->z,&T1,&T2);
	poly_mod_sqr(&P3->x,&T2);
	poly_mod_sqr(&T1,&T1);
	poly_mod_mult(&T2,&T1,&GRP->b);
	poly_add(&P3->x,&P3->x,&T2);
//	poly_mod_sqr(&T1,&P1->y);
//	poly_add(&T1,&T1,&P3->z);//if a==1 (i think it does)
//	poly_add(&T1,&T1,&T2);
//	poly_mod_mult(&P3->y,&P3->x,&T1);
//	poly_mod_mult(&T1,&T2,&P3->z);
//	poly_add(&P3->y,&P3->y,&T1);
}

void spec_point_LDadd(PECGROUP GRP, PECPOINT P3, PECPOINT P1, PECPOINT P2) {
	FIELD2M T1,T2,T3;

	if(spec_point_cmpeq(P2,&spec_point_Ainf)) {
		*P3=*P1; return;
	}
	if(spec_ppoint_cmpeq(P1,&spec_point_LDinf)) {
		*P3=*P2; P3->z=poly_one;
		return;
	}
	
	poly_mod_mult(&T1,&P1->z,&P2->x);
	poly_mod_sqr(&T2,&P1->z);
	poly_add(&P3->x,&P1->x,&T1);
	poly_mod_mult(&T1,&P1->z,&P3->x);
	poly_mod_mult(&T3,&T2,&P2->y);
	poly_add(&P3->y,&P1->y,&T3);

	if(poly_cmpeq(&P3->x,&poly_zero)) {
		if(poly_cmpeq(&P3->y,&poly_zero)) {
			spec_point_LDdbl(GRP,P3,P2);
		}else{
			*P3=spec_point_LDinf;
		}
		return;
	}

	poly_mod_sqr(&P3->z,&T1); 
	poly_mod_mult(&T3,&T1,&P3->y);
	poly_add(&T1,&T1,&T2);//if a==1 (which i think it does)
	poly_mod_sqr(&T2,&P3->x);
	poly_mod_mult(&P3->x,&T2,&T1);
	poly_mod_sqr(&T2,&P3->y);
	poly_add(&P3->x,&P3->x,&T2);
	poly_add(&P3->x,&P3->x,&T3);
	poly_mod_mult(&T2,&P2->x,&P3->z);
	poly_add(&T2,&T2,&P3->x);
	poly_mod_sqr(&T1,&P3->z);
	poly_add(&T3,&T3,&P3->z);
	poly_mod_mult(&P3->y,&T3,&T2);
	poly_add(&T2,&P2->x,&P2->y);
	poly_mod_mult(&T3,&T1,&T2);
	poly_add(&P3->y,&P3->y,&T3);
}
void spec_point_MLDadd(PECGROUP GRP, PFIELD2M x, PECPOINT P1, PECPOINT P2) {
	FIELD2M T1,T2,T3;

	if(spec_ppoint_cmpeq(P2,&spec_point_LDinf)) {
		return; //*P3=*P1; return;
	}
	if(spec_ppoint_cmpeq(P1,&spec_point_LDinf)) {
		*P1=*P2; //P3->z=poly_one;
		return;
	}
	
	poly_mod_mult(&T2,&P2->z,&P1->x);
	poly_mod_mult(&T1,&P1->z,&P2->x);
	poly_add(&T2,&T1,&T2);
	poly_mod_sqr(&P1->z,&T2);
	poly_mod_mult(&T2,&T1,&T2);
	poly_mod_sqr(&T1,&T1);
	poly_mod_mult(&P1->x,x,&P1->z);
	poly_add(&P1->x,&P1->x,&T2);
	poly_add(&P1->x,&P1->x,&T1);
}

void spec_point_LDkP(PECGROUP GRP, PECPOINT Q, PFIELD2M k, PECPOINT P) {
	int i;
	int Lminus1;
	int Lminus2mod8;
	unsigned char kw,ki;
	int bi,ei;
	FIELD2M Px=P->x;

//	FIELD2M P1,P2;
	ECPOINT Pn[2];

	*Q=spec_point_LDinf;

	Lminus1 = poly_degree(k);
	Lminus2mod8 = ((Lminus1-0) & 7);

	ei = Lminus1/ELEMENTBITS;
	bi = 15-((Lminus1%ELEMENTBITS)/8);

	DBG({
		printf("k= ");
		poly_print(k); printf("\n");
		printf("L-1 == %d\n",Lminus1);
		printf("ei : %d\n",ei);
		printf("bi : %d\n",bi);
		printf("(L-1 mod 8) == %d\n",Lminus2mod8);
	})

	kw = spu_extract((vec_uchar16)k->v[ei],bi);
	DBG({printf("kw (pre shift) = 0x%02x\n",kw);})
	kw<<=(7-Lminus2mod8);

	DBG({printf("kw = 0x%02x\n",kw);})

#define spec_point_LDadddbl(b) \
		spec_point_LDdbl(GRP,Q,Q); \
		if(b) spec_point_LDadd(GRP, Q, Q, P); 

	// do L-2 mod 8 steps first, so we can get to a byte boundry in k, then do a byte at a time 
	for(i=0;i<(Lminus2mod8+1);++i) {
		ki=(kw&0x80) >> 7;
		DBG({printf("ki = %d\n",ki);})
		spec_point_LDadddbl(ki);
		kw<<=1;
	}

	bi+=1; //go to next byte (first full byte)
	ei-=bi>>4; //if bi is 16, then ei goes to the next element
	bi&=0x0f; //bi = bi mod 16;

	while(ei>=0) {
		kw=spu_extract((vector unsigned char)(k->v[ei]),bi);
		DBG({printf("ei = %d\nbi = %d\nkw = 0x%02x\n",ei,bi,kw);})
		spec_point_LDadddbl((kw>>7)&1);
		spec_point_LDadddbl((kw>>6)&1);
		spec_point_LDadddbl((kw>>5)&1);
		spec_point_LDadddbl((kw>>4)&1);
		spec_point_LDadddbl((kw>>3)&1);
		spec_point_LDadddbl((kw>>2)&1);
		spec_point_LDadddbl((kw>>1)&1);
		spec_point_LDadddbl((kw)&1);
		bi+=1; //go to next byte in vector
		ei-=bi>>4; //if bi is 16, then ei goes to the next element
		bi&=0x0f; //bi = bi mod 16;
	}
//	DBG(printf("LdkP in projected coords:\n");spec_point_print(Pn+1);)
	spec_point_LD2A(GRP, Q, Q);
}

void spec_point_Madd(PECGROUP GRP, PFIELD2M x, PECPOINT P1, PECPOINT P2) {
	FIELD2M T1;
	FIELD2M T2;

	T1=*x;
	poly_mod_mult(&P1->x,&P1->x,&P2->z);
	poly_mod_mult(&P1->z,&P1->z,&P2->x);
	poly_mod_mult(&T2,&P1->x,&P1->z);
	poly_add(&P1->z,&P1->z,&P1->x);
	poly_mod_sqr(&P1->z,&P1->z);
	poly_mod_mult(&P1->x,&P1->z,&T1);
	poly_add(&P1->x,&P1->x,&T2);
}

// extracted from 3.40
/*
void spec_point_Mdbl(PECGROUP GRP, PECPOINT P) {
	FIELD2M T1,T2,T3;

	poly_mod_sqr(&T1,&P->x); //T1 = x^2
	poly_mod_sqr(&P->x,&T1); //x = x^4

	poly_mod_sqr(&T2,&P->z); //T2=Z^2
	poly_mod_sqr(&T3,&T2); //T3=Z^4
	poly_mod_mult(&T3,&T3,&GRP->b);//T3 = b*Z^4

	poly_add(&P->x,&P->x,&T3);//x+=bZ^4 = x^4 + bz^4

	poly_mod_mult(&P->z,&T1,&T2);//z=x^2 * z^2 = T1 * T2
	
}
*/

// from L-D paper
void spec_point_Mdbl(PECGROUP GRP, PECPOINT P) {
	FIELD2M T1;
/*
	T1=GRP->c;
	poly_mod_sqr(&P->x,&P->x);
	poly_mod_sqr(&P->z,&P->z);
	poly_mod_mult(&T1,&P->z,&T1);
	poly_mod_mult(&P->z,&P->z,&P->x);
	poly_mod_sqr(&T1,&T1);
	poly_mod_sqr(&P->x,&P->x);
	poly_add(&P->x,&P->x,&T1);
*/
// converted to NOT USE c (until i calculate it later);
	T1=GRP->c;
	poly_mod_sqr(&P->x,&P->x);
	poly_mod_sqr(&P->z,&P->z);
	poly_mod_mult(&T1,&P->z,&T1);
	poly_mod_mult(&P->z,&P->z,&P->x);
	poly_mod_sqr(&T1,&T1);
	poly_mod_sqr(&P->x,&P->x);
	poly_add(&P->x,&P->x,&T1);
}

// x,y = P2.x,z 
void spec_point_Mxy(PECGROUP GRP, PECPOINT P, PECPOINT P1, PECPOINT P2) {
	FIELD2M T1,T2,T3,T4,T5;

	DBG({
		FIELD2M Zinv;
		printf("Mxy CALLED with:\n");
		printf("P1:\n");
		spec_point_print(P1);
		if(!poly_cmpeq(P1->z,&poly_zero)) {
		printf("X/Z = ");
		poly_inv_mod(&Zinv,&P1->z);
		poly_mod_mult(&Zinv,&Zinv,&P1->x);
		poly_print(&Zinv);
		}
		printf("P2:\n");
		spec_point_print(P2);
		if(!poly_cmpeq(P2->z,&poly_zero)) {
		printf("X/Z = ");
		poly_inv_mod(&Zinv,&P2->z);
		poly_mod_mult(&Zinv,&Zinv,&P2->x);
		poly_print(&Zinv);
		}
	})

	if(poly_cmpeq(&P1->z,&poly_zero)) { P2->x=poly_zero; P2->z=poly_zero; return;}
	if(poly_cmpeq(&P2->z,&poly_zero)) { P2->x=P->x; poly_add(&P2->z,&P->x,&P->y); return; }

	T1=P->x;
	T2=P->y;
	poly_mod_mult(&T3,&P1->z,&P2->z);
	poly_mod_mult(&P1->z,&P1->z,&T1);
	poly_add(&P1->z,&P1->z,&P1->x);
	poly_mod_mult(&P2->z,&P2->z,&T1);
	poly_mod_mult(&P1->x,&P2->z,&P1->x);
	poly_add(&P2->z,&P2->z,&P2->x);
	poly_mod_mult(&P2->z,&P2->z,&P1->z);
	poly_mod_sqr(&T4,&T1);
	poly_add(&T4,&T4,&T2);
	poly_mod_mult(&T4,&T4,&T3);
	poly_add(&T4,&T4,&P2->z);
	poly_mod_mult(&T3,&T3,&T1);
	poly_inv_mod(&T5,&T3); // use T5 below instead of T3
	poly_mod_mult(&T4,&T5,&T4);
	poly_mod_mult(&P2->x,&P1->x,&T5);
	poly_add(&P2->z,&P2->x,&T1);
	poly_mod_mult(&P2->z,&P2->z,&T4);
	poly_add(&P2->z,&P2->z,&T2);
}

void spec_point_MkP(PECGROUP GRP, PECPOINT Q, PFIELD2M k, PECPOINT P) {
	int i;
	int Lminus1;
	int Lminus2mod8;
	unsigned char kw,ki;
	int bi,ei;
	FIELD2M Px=P->x;

//	FIELD2M P1,P2;
	ECPOINT Pn[2];

	Q->z=poly_one;
	if(poly_cmpeq(k,&poly_zero) || poly_cmpeq(&P->x,&poly_zero)) {
		Q->x=poly_zero; Q->y=poly_zero;
		return;
	}
	Lminus1 = poly_degree(k);
	Lminus2mod8 = ((Lminus1-1) & 7);
	if(Lminus1==0) {
		// degree 0 implies that k==1 since k!=0 if it gets here
		if(Q!=P) *Q=*P;
		return;
	}
	Lminus1-=1;

	ei = Lminus1/ELEMENTBITS;
	bi = 15-((Lminus1%ELEMENTBITS)/8);

	DBG({
		printf("k= ");
		poly_print(k); printf("\n");
		printf("L-1 == %d\n",Lminus1);
		printf("ei : %d\n",ei);
		printf("bi : %d\n",bi);
		printf("(L-2 mod 8)+1 == %d\n",Lminus2mod8);
	})

//	P1.x=P->x; P1.z=poly_one;
	Pn[0].x=P->x; Pn[0].z=poly_one;
//	poly_mod_sqr(&P2.z,&P->x); // Z_2 = x^2
	poly_mod_sqr(&Pn[1].z,&P->x); // Z_2 = x^2
//	poly_mod_sqr(&P2.x,&P2.z); // X_2 = Z_2^2 = x^4 (+ b -- next line)
	poly_mod_sqr(&Pn[1].x,&Pn[1].z); // X_2 = Z_2^2 = x^4 (+ b -- next line)
//	poly_add(&P2.x,GRP->b);
	poly_add(&Pn[1].x,&Pn[1].x,&GRP->b);
//	Pn[0]=*P; //DEBUG LD only REMOVE
//	spec_point_LDdbl(GRP, &Pn[1],&Pn[0]);//DEBUG LD only REMOVE
//	Pn[0]=*P; //DEBUG LD only REMOVE
//	spec_point_MLDdbl(GRP, &Pn[1],&Pn[0]);//DEBUG LD only REMOVE

	kw = spu_extract((vec_uchar16)k->v[ei],bi);
	kw<<=(7-Lminus2mod8);

	DBG({printf("kw = 0x%02x\n",kw);})

	// do L-2 mod 8 steps first, so we can get to a byte boundry in k, then do a byte at a time 
	for(i=0;i<Lminus2mod8+1;++i) {
		ki=(kw&0x80) >> 7;
		DBG({printf("ki = %d\n",ki);})
//		spec_point_Madd(GRP, &Px, Pn+(ki^1), Pn+ki);
//		spec_point_Mdbl(GRP, Pn+ki);
//
//		spec_point_LD2A(GRP,Pn+ki,Pn+ki);
//		spec_point_LDadd(GRP, Pn+(ki^1), Pn+(ki^1), Pn+ki);
//		spec_point_LDdbl(GRP, Pn+ki, Pn+ki);

		spec_point_MLDadd(GRP, &Px, Pn+(ki^1), Pn+ki);
		spec_point_MLDdbl(GRP, Pn+ki, Pn+ki);

		kw<<=1;
	}

	bi+=1; //no to next byte (first full byte)
	ei-=bi>>4; //if bi is 16, then ei goes to the next element
	bi&=0x0f; //bi = bi mod 16;

	while(ei>=0) {
		kw=spu_extract((vector unsigned char)(k->v[ei]),bi);
		DBG({printf("ei = %d\nbi = %d\nkw = 0x%02x\n",ei,bi,kw);})
#define spec_point_Madddbl(b) \
		spec_point_MLDadd(GRP, &Px, Pn+((b)^1), Pn+(b)); \
		spec_point_MLDdbl(GRP, Pn+(b), Pn+(b));
/*
#define spec_point_Madddbl(b) \
		spec_point_LDadd(GRP,Pn+((b)^1), Pn+((b)^1), Pn+(b)); \
		spec_point_LDdbl(GRP,Pn+(b),Pn+(b))
*/
		spec_point_Madddbl((kw>>7)&1);
		spec_point_Madddbl((kw>>6)&1);
		spec_point_Madddbl((kw>>5)&1);
		spec_point_Madddbl((kw>>4)&1);
		spec_point_Madddbl((kw>>3)&1);
		spec_point_Madddbl((kw>>2)&1);
		spec_point_Madddbl((kw>>1)&1);
		spec_point_Madddbl((kw)&1);
		bi+=1; //go to next byte in vector
		ei-=bi>>4; //if bi is 16, then ei goes to the next element
		bi&=0x0f; //bi = bi mod 16;
	}
//	DBG(printf("MkP in projected coords:\n");spec_point_print(Pn+1);)
	spec_point_Mxy(GRP, P, Pn, Pn+1);
	Q->x=Pn[1].x; Q->y=Pn[1].z; Q->z=poly_one;
//	spec_point_LD2A(GRP,Q,Pn);//LD ONLY REMOVE
}

void spec_point_MkP_be(PECGROUP GRP, PECPOINT Q, PFIELD2M k, PECPOINT P) {
	int i;
	int Lminus1;
	int Lminus2mod8;
	unsigned char kw,ki;
	int bi,ei;
	FIELD2M Px=P->x;

//	FIELD2M P1,P2;
	ECPOINT Pn[2];

	Q->z=poly_one;
	if(poly_cmpeq(k,&poly_zero) || poly_cmpeq(&P->x,&poly_zero)) {
		Q->x=poly_zero; Q->y=poly_zero;
		return;
	}
	Lminus1 = poly_degree_be(k);
	Lminus2mod8 = ((Lminus1-1) & 7);
	if(Lminus1==0) {
		// degree 0 implies that k==1 since k!=0 if it gets here
		if(Q!=P) *Q=*P;
		return;
	}
	Lminus1-=1;

	ei = MAXWORDS - (Lminus1/ELEMENTBITS) - 1;
	bi = 15-((Lminus1%ELEMENTBITS)/8);

	DBG({
		printf("k= ");
		poly_print(k); printf("\n");
		printf("L-1 == %d\n",Lminus1);
		printf("ei : %d\n",ei);
		printf("bi : %d\n",bi);
		printf("(L-2 mod 8)+1 == %d\n",Lminus2mod8);
	})

//	P1.x=P->x; P1.z=poly_one;
	Pn[0].x=P->x; Pn[0].z=poly_one;
//	poly_mod_sqr(&P2.z,&P->x); // Z_2 = x^2
	poly_mod_sqr(&Pn[1].z,&P->x); // Z_2 = x^2
//	poly_mod_sqr(&P2.x,&P2.z); // X_2 = Z_2^2 = x^4 (+ b -- next line)
	poly_mod_sqr(&Pn[1].x,&Pn[1].z); // X_2 = Z_2^2 = x^4 (+ b -- next line)
//	poly_add(&P2.x,GRP->b);
	poly_add(&Pn[1].x,&Pn[1].x,&GRP->b);
//	Pn[0]=*P; //DEBUG LD only REMOVE
//	spec_point_LDdbl(GRP, &Pn[1],&Pn[0]);//DEBUG LD only REMOVE
//	Pn[0]=*P; //DEBUG LD only REMOVE
//	spec_point_MLDdbl(GRP, &Pn[1],&Pn[0]);//DEBUG LD only REMOVE

	kw = spu_extract((vec_uchar16)k->v[ei],bi);
	kw<<=(7-Lminus2mod8);

	DBG({printf("kw = 0x%02x\n",kw);})

	// do L-2 mod 8 steps first, so we can get to a byte boundry in k, then do a byte at a time 
	for(i=0;i<Lminus2mod8+1;++i) {
		ki=(kw&0x80) >> 7;
		DBG({printf("ki = %d\n",ki);})
//		spec_point_Madd(GRP, &Px, Pn+(ki^1), Pn+ki);
//		spec_point_Mdbl(GRP, Pn+ki);
//
//		spec_point_LD2A(GRP,Pn+ki,Pn+ki);
//		spec_point_LDadd(GRP, Pn+(ki^1), Pn+(ki^1), Pn+ki);
//		spec_point_LDdbl(GRP, Pn+ki, Pn+ki);

		spec_point_MLDadd(GRP, &Px, Pn+(ki^1), Pn+ki);
		spec_point_MLDdbl(GRP, Pn+ki, Pn+ki);

		kw<<=1;
	}

	bi+=1; //no to next byte (first full byte)
	ei+=bi>>4; //if bi is 16, then ei goes to the next element
	bi&=0x0f; //bi = bi mod 16;

	while(ei<(signed int)MAXWORDS) {
		kw=spu_extract((vector unsigned char)(k->v[ei]),bi);
		DBG({printf("ei = %d\nbi = %d\nkw = 0x%02x\n",ei,bi,kw);})
/* //defined in non _be version 
#define spec_point_Madddbl(b) \
		spec_point_MLDadd(GRP, &Px, Pn+((b)^1), Pn+(b)); \
		spec_point_MLDdbl(GRP, Pn+(b), Pn+(b));
*/
/*
#define spec_point_Madddbl(b) \
		spec_point_LDadd(GRP,Pn+((b)^1), Pn+((b)^1), Pn+(b)); \
		spec_point_LDdbl(GRP,Pn+(b),Pn+(b))
*/
		spec_point_Madddbl((kw>>7)&1);
		spec_point_Madddbl((kw>>6)&1);
		spec_point_Madddbl((kw>>5)&1);
		spec_point_Madddbl((kw>>4)&1);
		spec_point_Madddbl((kw>>3)&1);
		spec_point_Madddbl((kw>>2)&1);
		spec_point_Madddbl((kw>>1)&1);
		spec_point_Madddbl((kw)&1);
		bi+=1; //go to next byte in vector
		ei+=bi>>4; //if bi is 16, then ei goes to the next element
		bi&=0x0f; //bi = bi mod 16;
	}
//	DBG(printf("MkP in projected coords:\n");spec_point_print(Pn+1);)
	spec_point_Mxy(GRP, P, Pn, Pn+1);
	Q->x=Pn[1].x; Q->y=Pn[1].z; Q->z=poly_one;
//	spec_point_LD2A(GRP,Q,Pn);//LD ONLY REMOVE
}

void spec_point_simulMkP_be(PECGROUP GRP, PECPOINT R, PFIELD2M k, PECPOINT P, PFIELD2M l, PECPOINT Q) {
	int i;
	int Lminus1,degL,degK;
	int Lminus2mod8;
	unsigned char kw,lw,ki;
	int bi,ei;
//	ECPOINT tbl[16];
	ECPOINT tbl[4][4];
//	PECPOINT tbl[2][2];
	ECPOINT PplusQ;
	FIELD2M zinv;

	*R=spec_point_LDinf;

	DBG({
		printf("Simul Mult:\nR:\n");
		spec_point_print(R);
		printf("Simul Mult:\nP:\n");
		spec_point_print(P);
		printf("Simul Mult:\nQ:\n");
		spec_point_print(Q);
	})

	/*calculate precompute table of tbl[i] = aP + bQ for a,b elementof {0,1,2,3} w=2,  i=4a+b*/
#if 1 
	tbl[0][0] = spec_point_Ainf;
	spec_point_LDadd(GRP,&tbl[0][1],&spec_point_LDinf,Q);
	spec_point_LDdbl(GRP,&tbl[0][2],&tbl[0][1]);
	spec_point_LDadd(GRP,&tbl[0][3],&tbl[0][2],Q);

	spec_point_LDadd(GRP,&tbl[1][0],&spec_point_LDinf,P);
	spec_point_LDdbl(GRP,&tbl[2][0],&tbl[1][0]);
	spec_point_LDadd(GRP,&tbl[3][0],&tbl[2][0],P);

	spec_point_LDadd(GRP,&tbl[1][1],&tbl[1][0],Q);
	spec_point_LDadd(GRP,&tbl[2][1],&tbl[2][0],Q);
	spec_point_LDadd(GRP,&tbl[3][1],&tbl[3][0],Q);

	spec_point_LDadd(GRP,&tbl[1][2],&tbl[1][1],Q);
	spec_point_LDdbl(GRP,&tbl[2][2],&tbl[1][1]);
	spec_point_LDadd(GRP,&tbl[3][2],&tbl[3][1],Q);

	spec_point_LDadd(GRP,&tbl[1][3],&tbl[1][2],Q);
	spec_point_LDadd(GRP,&tbl[2][3],&tbl[2][2],Q);
	spec_point_LDadd(GRP,&tbl[3][3],&tbl[3][2],Q);

	tbl[0][1] = *Q;
	spec_point_LD2A(GRP,&tbl[0][2],&tbl[0][2]);
	spec_point_LD2A(GRP,&tbl[0][3],&tbl[0][3]);
	tbl[1][0] = *P;
	spec_point_LD2A(GRP,&tbl[1][1],&tbl[1][1]);
	spec_point_LD2A(GRP,&tbl[1][2],&tbl[1][2]);
	spec_point_LD2A(GRP,&tbl[1][3],&tbl[1][3]);
	spec_point_LD2A(GRP,&tbl[2][0],&tbl[2][0]);
	spec_point_LD2A(GRP,&tbl[2][1],&tbl[2][1]);
	spec_point_LD2A(GRP,&tbl[2][2],&tbl[2][2]);
	spec_point_LD2A(GRP,&tbl[2][3],&tbl[2][3]);
	spec_point_LD2A(GRP,&tbl[3][0],&tbl[3][0]);
	spec_point_LD2A(GRP,&tbl[3][1],&tbl[3][1]);
	spec_point_LD2A(GRP,&tbl[3][2],&tbl[3][2]);
	spec_point_LD2A(GRP,&tbl[3][3],&tbl[3][3]);
#else
	tbl[0][0]=&spec_point_Ainf;
	tbl[1][0]=P;
	tbl[0][1]=Q;
	tbl[1][1]=&PplusQ;
	P->z=poly_one;
	spec_point_LDadd(GRP,&PplusQ,P,Q);
	spec_point_LD2A(GRP,&PplusQ,&PplusQ);
#endif

	DBG({
		ECPOINT T;
		spec_point_LD2A(GRP,&T,&tbl[0][3]);
		printf("Simul Mult:\nk: ");
		poly_print(k);
		printf("\nSimul Mult:\nl: ");
		poly_print(l);
		printf("\nSimul Mult:\ntbl[2][3]: ");
		spec_point_print(&tbl[0][3]);
		printf("\nSimul Mult:\ntbl[2][3]: ");
		spec_point_print(&T);
		printf("\n");
	})	


/*
	if(poly_cmpeq(k,&poly_zero) || poly_cmpeq(&P->x,&poly_zero)) {
		R->x=poly_zero; R->y=poly_zero;
		return;
	}
*/
	degL = poly_degree_be(l);
	degK = poly_degree_be(k);
	Lminus1 = degL>degK?degL:degK;
	Lminus1-=1;

	ei = MAXWORDS - (Lminus1/ELEMENTBITS) - 1;
	bi = 15-((Lminus1%ELEMENTBITS)/8);

	//start at this byte//bi+=1; //no to next byte (first full byte)
	ei+=bi>>4; //if bi is 16, then ei goes to the next element
	bi&=0x0f; //bi = bi mod 16;

	DBG(printf("Simul Mult: ei=%d, bi=%d\n",ei,bi);)

	while(ei<(signed int)MAXWORDS) {
		kw=spu_extract((vector unsigned char)(k->v[ei]),bi);
		lw=spu_extract((vector unsigned char)(l->v[ei]),bi);
		DBG({printf("ei = %d\nbi = %d\nkw = 0x%02x\n",ei,bi,kw);})
#if 1
#define spec_point_LDsimuladd(a,b) \
		DBG(printf("Accumulate a=%01x b=%01x\n",a,b);) \
		spec_point_LDdbl(GRP,R,R); \
		spec_point_LDdbl(GRP,R,R); \
		spec_point_LDadd(GRP,R,R,&tbl[a][b]);
#define spec_point_LDsimuladd_ms(s) \
		spec_point_LDsimuladd((kw>>(2*(s-1)))&0x3,(lw>>(2*(s-1)))&0x3);
		spec_point_LDsimuladd_ms(4);	
		spec_point_LDsimuladd_ms(3);	
		spec_point_LDsimuladd_ms(2);	
		spec_point_LDsimuladd_ms(1);	
#else
#define spec_point_LDsimuladd(a,b) \
		DBG(printf("Accumulate a=%01x b=%01x\n",a,b);) \
		spec_point_LDdbl(GRP,R,R); \
		spec_point_LDadd(GRP,R,R,tbl[a][b]);
#define spec_point_LDsimuladd_ms(s) \
		spec_point_LDsimuladd((kw>>((s-1)))&0x1,(lw>>((s-1)))&0x1);
	
		spec_point_LDsimuladd_ms(8);	
		spec_point_LDsimuladd_ms(7);	
		spec_point_LDsimuladd_ms(6);	
		spec_point_LDsimuladd_ms(5);	
		spec_point_LDsimuladd_ms(4);	
		spec_point_LDsimuladd_ms(3);	
		spec_point_LDsimuladd_ms(2);	
		spec_point_LDsimuladd_ms(1);	
#endif
		bi+=1; //go to next byte in vector
		ei+=bi>>4; //if bi is 16, then ei goes to the next element
		bi&=0x0f; //bi = bi mod 16;
	}
//	DBG(printf("MkP in projected coords:\n");spec_point_print(Pn+1);)
//	spec_point_Mxy(GRP, P, Pn, Pn+1);
//	Q->x=Pn[1].x; Q->y=Pn[1].z; Q->z=poly_one;
//	spec_point_LD2A(GRP,Q,Pn);//LD ONLY REMOVE

	spec_point_LD2A(GRP,R,R);
//	poly_inv_mod(&zinv,&R->z);
//	poly_mod_mult(&R->x,&R->x,&zinv);
//	poly_mod_mult(&R->y,&R->y,&zinv);
//	R->z=poly_one;
}
void spec_point_simulMkP_be_w4small(PECGROUP GRP, PECPOINT R, PFIELD2M k, PECPOINT P, PFIELD2M l, PECPOINT Q) {
	int i;
	int Lminus1,degL,degK;
	int Lminus2mod8;
	unsigned char kw,lw,ki;
	int bi,ei;
//	ECPOINT tbl[16];
	ECPOINT tbl[2][16];
//	PECPOINT tbl[2][2];
	ECPOINT PplusQ;
	FIELD2M zinv;

	*R=spec_point_LDinf;

	DBG({
		printf("Simul Mult:\nR:\n");
		spec_point_print(R);
		printf("Simul Mult:\nP:\n");
		spec_point_print(P);
		printf("Simul Mult:\nQ:\n");
		spec_point_print(Q);
	})

	/*calculate precompute table of tbl[i] = aP + bQ for a,b elementof {0,1,2,3} w=2,  i=4a+b*/
#if 1 
	tbl[0][0] = tbl[1][0] = spec_point_Ainf;
	spec_point_LDadd(GRP,&tbl[0][1],&spec_point_LDinf,P);
	for(i=2;i<16;i+=2) {
		spec_point_LDdbl(GRP,&tbl[0][i],&tbl[0][i/2]);
		spec_point_LDadd(GRP,&tbl[0][i+1],&tbl[0][i],P);
	}
	tbl[0][1] = *P;
	for(i=2;i<16;++i) {
		spec_point_LD2A(GRP,&tbl[0][i],&tbl[0][i]);
	}
	spec_point_LDadd(GRP,&tbl[1][1],&spec_point_LDinf,Q);
	for(i=2;i<16;i+=2) {
		spec_point_LDdbl(GRP,&tbl[1][i],&tbl[1][i/2]);
		spec_point_LDadd(GRP,&tbl[1][i+1],&tbl[1][i],Q);
	}
	tbl[1][1] = *Q;
	for(i=2;i<16;++i) {
		spec_point_LD2A(GRP,&tbl[1][i],&tbl[1][i]);
	}

#else
	tbl[0][0]=&spec_point_Ainf;
	tbl[1][0]=P;
	tbl[0][1]=Q;
	tbl[1][1]=&PplusQ;
	P->z=poly_one;
	spec_point_LDadd(GRP,&PplusQ,P,Q);
	spec_point_LD2A(GRP,&PplusQ,&PplusQ);
#endif

	DBG({
		ECPOINT T;
		spec_point_LD2A(GRP,&T,&tbl[0][3]);
		printf("Simul Mult:\nk: ");
		poly_print(k);
		printf("\nSimul Mult:\nl: ");
		poly_print(l);
		printf("\nSimul Mult:\ntbl[2][3]: ");
		spec_point_print(&tbl[0][3]);
		printf("\nSimul Mult:\ntbl[2][3]: ");
		spec_point_print(&T);
		printf("\n");
	})	


/*
	if(poly_cmpeq(k,&poly_zero) || poly_cmpeq(&P->x,&poly_zero)) {
		R->x=poly_zero; R->y=poly_zero;
		return;
	}
*/
	degL = poly_degree_be(l);
	degK = poly_degree_be(k);
	Lminus1 = degL>degK?degL:degK;
	Lminus1-=1;

	ei = MAXWORDS - (Lminus1/ELEMENTBITS) - 1;
	bi = 15-((Lminus1%ELEMENTBITS)/8);

	//start at this byte//bi+=1; //no to next byte (first full byte)
	ei+=bi>>4; //if bi is 16, then ei goes to the next element
	bi&=0x0f; //bi = bi mod 16;

	DBG(printf("Simul Mult: ei=%d, bi=%d\n",ei,bi);)

	while(ei<(signed int)MAXWORDS) {
		kw=spu_extract((vector unsigned char)(k->v[ei]),bi);
		lw=spu_extract((vector unsigned char)(l->v[ei]),bi);
		DBG({printf("ei = %d\nbi = %d\nkw = 0x%02x\n",ei,bi,kw);})
#if 1
#define spec_point_LDsimuladd_big(a,b) \
		DBG(printf("Accumulate a=%01x b=%01x\n",a,b);) \
		spec_point_LDdbl(GRP,R,R); \
		spec_point_LDdbl(GRP,R,R); \
		spec_point_LDdbl(GRP,R,R); \
		spec_point_LDdbl(GRP,R,R); \
		spec_point_LDadd(GRP,R,R,&tbl[0][a]); \
		spec_point_LDadd(GRP,R,R,&tbl[1][b]);
#define spec_point_LDsimuladd_ms_big(s) \
		spec_point_LDsimuladd_big((kw>>(4*(s-1)))&0xF,(lw>>(4*(s-1)))&0xF);
		spec_point_LDsimuladd_ms_big(2);	
		spec_point_LDsimuladd_ms_big(1);	
#else
#define spec_point_LDsimuladd(a,b) \
		DBG(printf("Accumulate a=%01x b=%01x\n",a,b);) \
		spec_point_LDdbl(GRP,R,R); \
		spec_point_LDadd(GRP,R,R,tbl[a][b]);
#define spec_point_LDsimuladd_ms(s) \
		spec_point_LDsimuladd((kw>>((s-1)))&0x1,(lw>>((s-1)))&0x1);
	
		spec_point_LDsimuladd_ms(8);	
		spec_point_LDsimuladd_ms(7);	
		spec_point_LDsimuladd_ms(6);	
		spec_point_LDsimuladd_ms(5);	
		spec_point_LDsimuladd_ms(4);	
		spec_point_LDsimuladd_ms(3);	
		spec_point_LDsimuladd_ms(2);	
		spec_point_LDsimuladd_ms(1);	
#endif
		bi+=1; //go to next byte in vector
		ei+=bi>>4; //if bi is 16, then ei goes to the next element
		bi&=0x0f; //bi = bi mod 16;
	}
//	DBG(printf("MkP in projected coords:\n");spec_point_print(Pn+1);)
//	spec_point_Mxy(GRP, P, Pn, Pn+1);
//	Q->x=Pn[1].x; Q->y=Pn[1].z; Q->z=poly_one;
//	spec_point_LD2A(GRP,Q,Pn);//LD ONLY REMOVE

	spec_point_LD2A(GRP,R,R);
//	poly_inv_mod(&zinv,&R->z);
//	poly_mod_mult(&R->x,&R->x,&zinv);
//	poly_mod_mult(&R->y,&R->y,&zinv);
//	R->z=poly_one;
}

#define call_mpm_mod(r,x,n) \
	mpm_mod((vector unsigned int*)(r).v,(vector unsigned int*)(x).v,MAXWORDS,(vector unsigned int*)(n).v,MAXWORDS)
#define call_mpm_mod_mul(r,x,y,n) \
	mpm_mul((vector unsigned int*)(_t.v),(vector unsigned int*)(x).v,MAXWORDS,(vector unsigned int*)(y).v,MAXWORDS); \
	mpm_mod((vector unsigned int*)(r).v,(vector unsigned int*)(_t.v),2*MAXWORDS,(vector unsigned int*)(n).v,MAXWORDS)
//this works because the overflow to the 234th bit is still under 256...
#define call_mpm_mod_add(r,x,y,n) \
	mpm_add((vector unsigned int*)(r).v,(vector unsigned int*)(x).v,(vector unsigned int*)(y).v,MAXWORDS); \
	call_mpm_mod(r,r,n)
#define call_mpm_mod_mul_inv(r,x,n) \
	mpm_mul_inv((vector unsigned int*)(r).v,(vector unsigned int*)(n).v,(vector unsigned int*)(x).v,MAXWORDS)
#define call_mpm_mod_div(r,x,y,n) \
	call_mpm_mod_mul_inv(_t1,y,n); \
	call_mpm_mod_mul(r,x,_t1,n)
#define call_mpm_cmpeq(a,b) \
	mpm_cmpeq((vector unsigned int*)(a).v,(vector unsigned int*)(b).v,MAXWORDS)


int spec_ecdsa_verify(PECGROUP GRP, PECPOINT Q, PFIELD2M e, PFIELD2M r, PFIELD2M s) {
	FIELD2M w,u1,u2,v;
	ECPOINT T,X;
	DBLFIELD2M _t,_t1;

	DBG({
		printf("ecdsa_verify:\n");
		printf("Q:\n"); spec_point_print(Q);
		printf("e:\n"); poly_print(e); printf("\n");
		printf("r:\n"); poly_print(r); printf("\n");
		printf("s:\n"); poly_print(s); printf("\n");
	})

	//verify r,s range //return -3
	
	call_mpm_mod_mul_inv(w,*s,GRP->nbe);
	call_mpm_mod_mul(u1,*e,w,GRP->nbe);
	call_mpm_mod_mul(u2,*r,w,GRP->nbe);
#if 0
	//spec_point_MkP(GRP,&X,&u1,&GRP->G);
	//spec_point_MkP(GRP,&T,&u2,Q);
	spec_point_MkP_be(GRP,&X,&u1,&GRP->G);
	spec_point_MkP_be(GRP,&T,&u2,Q);
	X.z=poly_one;
	spec_point_LDadd(GRP,&X,&X,&T);
	spec_point_LD2A(GRP, &X, &X);

	if(spec_point_cmpeq(&X,&spec_point_Ainf)) {
		return(-2);
	}
#else
	spec_point_simulMkP_be(GRP,&X,&u1,&GRP->G,&u2,Q);
#endif

	call_mpm_mod(v,X.x,GRP->nbe);
	
	if(call_mpm_cmpeq(v,*r)) {
		return(1);
	}else{
		return(0);
	}

}

int spec_ecdsa_sign(PECGROUP GRP, PFIELD2M kbe, PFIELD2M d, PFIELD2M e, PFIELD2M r, PFIELD2M s) {
	ECPOINT kP;
	DBLFIELD2M _t,_t1,t,kinv;

	DBG({
		FIELD2M a=field2m_init_ull4be(0,0,0,B);
		FIELD2M b=field2m_init_ull4be(0,0,0,7);
		FIELD2M c;
		FIELD2M cc;
		printf("ecdsa_sign:\n");
		printf("simple mod tests first\n");
		call_mpm_mod(c,a,b);
		printf("c = a % b ==? 4:\n"); poly_print(&c); printf("\n");
		a=c;
		cc=c;
		call_mpm_mod_mul(c,a,cc,b);
		printf("c = c^2 % b ==? 2:\n"); poly_print(&c); printf("\n");
		
		printf("kbe:\n"); poly_print(kbe); printf("\n");
		printf("d:\n"); poly_print(d); printf("\n");
		printf("e:\n"); poly_print(e); printf("\n");
	});

	call_mpm_mod(*kbe,*kbe,GRP->nbe);
	if(call_mpm_cmpeq(*kbe,poly_zero)) {
		return(-2); //bad choice for k
	}

	//spec_point_MkP(GRP,&kP,kle,&GRP->G);
	spec_point_MkP_be(GRP,&kP,kbe,&GRP->G);
	call_mpm_mod(*r,kP.x,GRP->nbe);
//	mpm_mod((vector unsigned int*)r->v,(vector unsigned int*)kP.x.v,MAXWORDS,(vector unsigned int*)GRP->n.v,MAXWORDS);
//static __inline void _mpm_mod(vector unsigned int *r, const vector unsigned int *a, int asize, const vector unsigned int *b, int bsize)

	if(call_mpm_cmpeq(*r,poly_zero)) {
		return(-1); //bad choice for k
	}


	call_mpm_mod_mul(t,*r,*d,GRP->nbe);
	call_mpm_mod_add(t,t,*e,GRP->nbe);
	call_mpm_mod_div(*s,t,*kbe,GRP->nbe);

	DBG({
		printf("r:\n"); poly_print(r); printf("\n");
		printf("s:\n"); poly_print(s); printf("\n");
	})
	return(0);
}


